# Replication file 03b_Case_level_AUC_and_table_creation
# Article: Counterfactual Coercion: Could harsher sanctions against Russia have prevented the worst?
# Authors: Thies Niemeier, Gerald Schneider

# set appropriate working directory
# setwd()

# install packages
# install.packages("readxl")
# install.packages("tidyverse")
# install.packages("pROC")

library(readxl)
library(tidyverse)
library(pROC)

# Success AUC computations
rm(list=ls())

## EU ####
EU <- read_excel("Supplemental_Material/EU_Case_level_Success_Results.xls", 
                 sheet = "Sheet1")
EU <- EU %>% filter(!is.na(sucEU))

# AUC among observed cases
EU <- EU %>%
  mutate(ordsuccess_test = if_else(ordsuccess_test==0, "X0",
                                   if_else(ordsuccess_test==1, "X1", "X2")))

response_sucEU <- as.factor(EU$ordsuccess_test)
predictor_sucEU <- as.matrix(data.frame(
  "0" = EU$sucEU1,
  "1" = EU$sucEU2,
  "2" = EU$sucEU3
))

EU_suc <- multiclass.roc(response_sucEU, predictor_sucEU)
# AUC for success of EU sanctions across 3 classes: 0.69
EU_suc

# remove irrelevant objects
rm(EU_suc, predictor_sucEU, response_sucEU)

## US ####
US <- read_excel("Supplemental_Material/US_Case_level_Success_Results.xls")
US <- US %>% filter(!is.na(sucUS))

# AUC (ROC for ordinal variables) among observed cases
US <- US %>%
  mutate(ordsuccess_test = if_else(ordsuccess_test==0, "X0",
                                   if_else(ordsuccess_test==1, "X1", "X2")))

response_sucUS <- as.factor(US$ordsuccess_test)
predictor_sucUS <- as.matrix(data.frame(
  "0" = US$sucUS1,
  "1" = US$sucUS2,
  "2" = US$sucUS3
))

US_suc <- multiclass.roc(response_sucUS, predictor_sucUS)
# AUC for success of US sanctions across 3 classes: .73
US_suc

# remove irrelevant objects
rm(US_suc, predictor_sucUS, response_sucUS)

# Create tables for main text, supplemental material
# Supplemental material
EU <- EU %>%
  select(sender, year, cname, caseid, HSEscore, sucEU, sucEUlow, sucEUmed, sucEUhigh, sucEUhigher) %>%
  mutate(caseid = as.character(caseid),
         sucEU = round(sucEU, digits = 2),
         sucEUlow = round(sucEUlow, digits = 2),
         sucEUmed = round(sucEUmed, digits = 2),
         sucEUhigh = round(sucEUhigh, digits = 2),
         sucEUhigher = round(sucEUhigher, digits = 2))

US<- US %>%
  select(sender, year, cname, caseid, HSEscore, sucUS, sucUSlow, sucUSmed, sucUShigh, sucUShigher) %>%
  mutate(caseid = as.character(caseid),
         sucUS = round(sucUS, digits = 2),
         sucUSlow = round(sucUSlow, digits = 2),
         sucUSmed = round(sucUSmed, digits = 2),
         sucUShigh = round(sucUShigh, digits = 2),
         sucUShigher = round(sucUShigher, digits = 2))

# Choice of cases
Case_choice <- EU %>%
  mutate(both_sanction = if_else(caseid %in% US$caseid, 1, 0)) %>%
  filter(both_sanction == 1) %>%
  mutate(chosen_cases = if_else(sucEU == min(sucEU), 1,
                                if_else(sucEU == max(sucEU), 1, 0)))
# Mean success prediction for EU: 0.799
mean(Case_choice$sucEU)
# Cases closest to 0.799 are Russia and Burundi
Case_choice %>%
  arrange(sucEU)

# Filter for chosen cases
Case_choice <- Case_choice %>%
  mutate(chosen_cases = if_else(caseid %in% c(2014022701, 2015050801), 1, chosen_cases)) %>%
  filter(chosen_cases == 1) 

# add US, rename variables
Case_choice <- Case_choice %>%
  left_join(US %>% select(caseid, sucUS, sucUSlow, sucUSmed, sucUShigh, sucUShigher), by = "caseid") %>%
  rename(Sender = sender,
         Year = year,
         Country_Name = cname,
         Case_ID = caseid,
         HSE = HSEscore,
         EU_Success_at_observed_int = sucEU,
         EU_Success_at_int_1 = sucEUlow,
         EU_Success_at_int_2 = sucEUmed,
         EU_Success_at_int_3 = sucEUhigh,
         EU_Success_at_int_4 = sucEUhigher,
         US_Success_at_observed_int = sucUS,
         US_Success_at_int_1 = sucUSlow,
         US_Success_at_int_2 = sucUSmed,
         US_Success_at_int_3 = sucUShigh,
         US_Success_at_int_4 = sucUShigher)

write.csv(Case_choice, file = "Main_Article/Counterfactual_Success_Predictions_4_Cases.csv")

# Give out table of success predictions for supplemental material
EU <- read_excel("Supplemental_Material/EU_Case_level_Success_Results.xls", 
                 sheet = "Sheet1")
US <- read_excel("Supplemental_Material/US_Case_level_Success_Results.xls")
colnames(US) <- colnames(EU)
Cases <- rbind(EU, US)

Cases <- Cases %>%
  select(sender, year, cname, caseid, HSEscore, sucEU, sucEUlow, sucEUmed, sucEUhigh, sucEUhigher) %>%
  mutate(caseid = as.character(caseid),
         sucEU = round(sucEU, digits = 2),
         sucEUlow = round(sucEUlow, digits = 2),
         sucEUmed = round(sucEUmed, digits = 2),
         sucEUhigh = round(sucEUhigh, digits = 2),
         sucEUhigher = round(sucEUhigher, digits = 2)) %>%
  rename(Sender = sender,
         Year = year,
         Country_Name = cname,
         Case_ID = caseid,
         HSE = HSEscore,
         Success_at_observed_int = sucEU,
         Success_at_int_1 = sucEUlow,
         Success_at_int_2 = sucEUmed,
         Success_at_int_3 = sucEUhigh,
         Success_at_int_4 = sucEUhigher)

# give out table of counterfactual success predictions across all cases
write.csv(Case_choice, file = "Supplemental_Material/Counterfactual_Success_Predictions_across_Cases.csv")

